# Gilardi, Fabrizio, "The Temporary Importance of Role Models for Women's Political Representation", American Journal of Political Science
# Code to replicate Figure 1 (Effect of an additional woman elected in other municipalities on the number of women candidates)
# gilardi@ipz.uzh.ch, 2014-06-23

# Set working directory
setwd("../Data/")

# Load packages
library(Zelig)

# Load data
d <- read.csv("dataset-full.csv")


# Function to compute mode
Mode <- function(x) {
 ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}


# Estimate models

out.2.z <- glm(n.women.cand ~ sl.autob.n.women.elected + n.women.elected.lagged + n.women.cand.previously.elected + n.women.cand.previously.elected + n.men.cand.previously.elected + n.women.cand.lagged + n.people.elected + n.cand + referenda + wom.kr + log.dist.zh + log.pop +  kr.svp + steuerfuss + agglo + land
	+ sl.autob.n.women.elected*y1978
	+ sl.autob.n.women.elected*y1982
	+ sl.autob.n.women.elected*y1986
	+ sl.autob.n.women.elected*y1990
	+ sl.autob.n.women.elected*y1994
	+ sl.autob.n.women.elected*y1998
	+ sl.autob.n.women.elected*y2002
	+ sl.autob.n.women.elected*y2006
	+ sl.autob.n.women.elected*y2010
	, data=d, family="poisson")
summary(out.2.z)


mm.2 <- model.frame(out.2.z)

z.out.2 <- zelig(n.women.cand ~ sl.autob.n.women.elected + n.women.elected.lagged + n.women.cand.previously.elected + n.women.cand.previously.elected + n.men.cand.previously.elected + n.women.cand.lagged + n.people.elected + n.cand + referenda + wom.kr + log.dist.zh + log.pop +  kr.svp + steuerfuss + agglo + land
	+ sl.autob.n.women.elected*y1978
	+ sl.autob.n.women.elected*y1982
	+ sl.autob.n.women.elected*y1986
	+ sl.autob.n.women.elected*y1990
	+ sl.autob.n.women.elected*y1994
	+ sl.autob.n.women.elected*y1998
	+ sl.autob.n.women.elected*y2002
	+ sl.autob.n.women.elected*y2006
	+ sl.autob.n.women.elected*y2010
	, model="poisson", data=mm.2)
summary(z.out.2)


# Vector of election years
yr <- seq(1974, 2010, 4)

# Significance level
sig <- qnorm(0.975)

# Empty matrix to store expected values
pr.l.2 <- pr.h.2 <- pr.d.2 <- matrix(NA, length(yr), 4)
pr.l.2[,1] <- pr.h.2[,1] <- pr.d.2[,1] <- yr

# Empty vetor to store baseline values
baseline.2 <- c()


# Compute expected values fot each election year

for(i in 1:length(yr)){

# Compute baseline value of the spatial lag (modal value)
sl.l <- round(Mode(d$sl.autob.n.women.elected[d$year==yr[i]]), 2)
baseline.2[i] <- sl.l

# Increase baseline by 0.05 (1 woman out of 20 neighbors)
sl.h <- sl.l + 0.05

# Vector of explanatory variables with baseline value for spatial lag
x.l <- setx(z.out.2, sl.autob.n.women.elected=sl.l, n.women.elected.lagged=round(Mode(d$n.women.elected.lagged[d$year==yr[i]]), 1), n.women.cand.lagged=round(Mode(d$n.women.cand.lagged[d$year==yr[i]]), 1), n.women.cand.previously.elected=round(Mode(d$n.women.cand.previously.elected[d$year==yr[i]]), 1), n.men.cand.previously.elected=round(Mode(d$n.men.cand.previously.elected[d$year==yr[i]]), 1),
	y1978 = unique(d$y1978[d$year==yr[i]]),
	y1982 = unique(d$y1982[d$year==yr[i]]),
	y1986 = unique(d$y1986[d$year==yr[i]]),
	y1990 = unique(d$y1990[d$year==yr[i]]),
	y1994 = unique(d$y1994[d$year==yr[i]]),
	y1998 = unique(d$y1998[d$year==yr[i]]),
	y2002 = unique(d$y2002[d$year==yr[i]]),
	y2006 = unique(d$y2006[d$year==yr[i]]),
	y2010 = unique(d$y2010[d$year==yr[i]])
	)

# Vector of explanatory variables with baseline value + 0.05 for spatial lag
x.h <- setx(z.out.2, sl.autob.n.women.elected=sl.h, n.women.elected.lagged=round(Mode(d$n.women.elected.lagged[d$year==yr[i]]), 1), n.women.cand.lagged=round(Mode(d$n.women.cand.lagged[d$year==yr[i]]), 1), n.women.cand.previously.elected=round(Mode(d$n.women.cand.previously.elected[d$year==yr[i]]), 1),  n.men.cand.previously.elected=round(Mode(d$n.men.cand.previously.elected[d$year==yr[i]]), 1),
	y1978 = unique(d$y1978[d$year==yr[i]]),
	y1982 = unique(d$y1982[d$year==yr[i]]),
	y1986 = unique(d$y1986[d$year==yr[i]]),
	y1990 = unique(d$y1990[d$year==yr[i]]),
	y1994 = unique(d$y1994[d$year==yr[i]]),
	y1998 = unique(d$y1998[d$year==yr[i]]),
	y2002 = unique(d$y2002[d$year==yr[i]]),
	y2006 = unique(d$y2006[d$year==yr[i]]),
	y2010 = unique(d$y2010[d$year==yr[i]])
	)

# Simulate values for baseline
s.out.l <- sim(z.out.2, x=x.l)
# Simulate values for baseline + 0.05
s.out.h <- sim(z.out.2, x=x.h)
# Simulate difference
s.out.d <- sim(z.out.2, x=x.l, x1=x.h)

# Expected value with baseline
pr.l.2[i,2] <- summary(s.out.l)$stats$"Expected Values: E(Y|X)"[1]
# Lower bound of confidence interval with baseline
pr.l.2[i,3] <- pr.l.2[i,2] - sig*summary(s.out.l)$stats$"Expected Values: E(Y|X)"[2]
# Upper bound of confidence interval with baseline
pr.l.2[i,4] <- pr.l.2[i,2] + sig*summary(s.out.l)$stats$"Expected Values: E(Y|X)"[2]

# Expected value with baseline + 0.05
pr.h.2[i,2] <- summary(s.out.h)$stats$"Expected Values: E(Y|X)"[1]
# Lower bound of confidence interval with baseline 1 0.05
pr.h.2[i,3] <- pr.h.2[i,2] - sig*summary(s.out.h)$stats$"Expected Values: E(Y|X)"[2]
# Upper bound of confidence interval with baseline + 0.05
pr.h.2[i,4] <- pr.h.2[i,2] + sig*summary(s.out.h)$stats$"Expected Values: E(Y|X)"[2]

# First difference
pr.d.2[i,2] <- summary(s.out.d)$stats$"First Differences: E(Y|X1) - E(Y|X)"[1]
# Lower bound of confidence interval for first difference
pr.d.2[i,3] <- pr.d.2[i,2] - sig*summary(s.out.d)$stats$"First Differences: E(Y|X1) - E(Y|X)"[2]
# Upper bound of confidence interval for first difference
pr.d.2[i,4] <- pr.d.2[i,2] + sig*summary(s.out.d)$stats$"First Differences: E(Y|X1) - E(Y|X)"[2]

}



#pdf(file="Figure-2.pdf", paper="special", width=7, height=5.5)
par(mar=c(3.5,3.5,2,1), mgp=c(2.5,0.8,0), cex.axis=1.2, cex.lab=1.2, font.main=1)
plot(yr, pr.d.2[,2], type="n", ylim=range(pr.d.2[,2:4]), xlab="", ylab="Expected increase in nr of women cand.", axes=F)
axis(1, at=yr)
axis(2)
abline(h=0, lty=2)
points(yr, pr.d.2[,2], type="b", cex=1.5)
segments(yr, pr.d.2[,3], yr, pr.d.2[,4])
legend("topright", pch=c(), legend="Effect of an additional woman elected\nin NEARBY municipalities (t-1)", pt.cex=1.5, cex=1.2, bty="n")
mtext(baseline.2*20, at=yr, side=1, padj=3, cex=0.9)
mtext("Baseline", at=1974, side=1, padj=3, adj=1.4, cex=0.9)
#dev.off()







